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ABSTRACT 

In the intracluster medium (ICM) of galaxy clusters, heat and momentum are trans- 
ported almost entirely along (but not across) magnetic field lines. We perform the first 
fully self-consistent Braginskii-MHD simulations of galaxy clusters including both of 
these effects. Specifically, we perform local and global simulations of the magnetother- 
mal instability (MTI) and the heat-flux-driven buoyancy instability (HBI) and assess 
the effects of viscosity on their saturation and astrophysical implications. We find that 
viscosity has only a modest effect on the saturation of the MTI. As in previous calcula- 
tions, we find that the MTI can generate nearly sonic turbulent velocities in the outer 
parts of galaxy clusters, although viscosity somewhat suppresses the magnetic field am- 
plification. At smaller radii in cool-core clusters, viscosity can decrease the linear growth 
rates of the HBI. However, it has less of an effect on the HBFs nonlinear saturation, in 
part because three-dimensional interchange motions (magnetic flux tubes slipping past 
each other) are not damped by anisotropic viscosity. In global simulations of cool core 
clusters, we show that the HBI robustly inhibits radial thermal conduction and thus pre- 
cipitates a cooling catastrophe. The effects of viscosity are, however, more important for 
higher entropy clusters. We argue that viscosity can contribute to the global transition 
of cluster cores from cool-core to non cool-core states: additional sources of intracluster 
turbulence, such as can be produced by AGN feedback or galactic wakes, suppress the 
HBI, heating the cluster core by thermal conduction; this makes the ICM more viscous, 
which slows the growth of the HBI, allowing further conductive heating of the cluster 
core and a transition to a non cool-core state. 

Key words: convection — galaxies: clusters: intracluster medium — instabilities — 
turbulence — X-rays: galaxies: clusters 



1 INTRODUCTION 

Clusters of galaxies are the largest gravitationally bound ob- 
jects in the universe, and as such, they potentially provide 
sensitive tests of cosmological parameters. They are filled with 
a hot, dilute, magnetized plasma, the intracluster medium 
(ICM), that emits copious X-rays. Observations of the ICM 
provide an interesting and unique window on problems rang- 
ing from constraining dark energy to understanding the ac- 
cretion and feedback processes for supermassive black holes. 

The plasma in the ICM has temperatures ranging from 
1-15 keV and number densities from 1CP 4 to 1CP 1 cm -3 . This 
dilute plasma has a mag netic field that has be en estimated 
to range from 0.1-10 fj,G (|Carilli fc Tavlorll2002h . With these 
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parameters the mean free path of electrons along the magnetic 
field line is > 10 12 times larger than the gyroradius at all 
radii in the ICM. Ions have a similar separation of scales. The 
mean free path is, however, always shorter than the local scale 
height; therefore, a fluid description of the plasma (as opposed 
to a collisionless description) is appropriate and th e ICM can 
be d escribed by the Braginskii-MHD equations l|Braginskiil 
1965). These equations are the standard ideal MHD equations 
supplemented with anisotropic conduction due to the electron 
heat flux and anisotropic momentum transport due to the 
ion viscosity along the magnetic field. In the ICM, transport 
perpendicular to the local magnetic field is negligible. 

As a result of the anisotropic heat transport in the 
ICM, the Schwarzschild criterion for convective instability — 
that the entropy increase in the direction of gravity — is re- 
placed by a criterion on temperature. In recent years, two 
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buoyant instabilities have been discovered that drive con- 
vection with the temperature gradient as a source of free 
energy. The first instabili ty, the magne tothermal instability 
(MTI), was described in iBalbusI (2000) and has been sim- 
ulated in two and three di mensions ( Parrish fc Stone] 120051 . 
120071 ; iMcCourt et al.ll2011al ). The MTI is unstable when the 
temperature gradient and gravity are in the same direction 
and grows fastest for a magnetic field perpendicular to grav- 
ity. The MTI operates in the outskirts of galaxy clusters and 
has been found to drive vigorous convection that can pro- 
vide over 30% of th e pressure support near the virial radius 
(|Parrish et al.l l201ll ). The second instability, the heat-flux- 
driven buoyancy instability (HBI) was described in lQuataertl 
(|2008l) and has been simulated in local simulations in 2D and 
3D i|Parrish fc Quataertll2008h . The HBI is unstable when the 
temperature gradient and gravity point in opposite directions 
and has the fastest growth for a magnetic field parallel to 
gravity. The HBI operates in the centers of cool-core clusters 
and saturates by reorienting the magnetic field to be perpen- 
dicular to gravity, greatly reducing the effecti ve radial con- 
ductivity and hastening a coo ling catastrophe (|Parrish et al.l 
120091 : iBogdanovic et all 120091 ). The interaction of the HBI 
with turbulence can help explain the bimod ality observed 
between cool core and non-c ool core clusters (Pa rrish et al.l 
|2010| ; iRuszkowski fc Ohll2010l ) . 

None of these previous numerical studies have self- 
consistently included viscosity because the ratio of viscous 
to thermal diffusion, the Prandtl numberQ for a hydrogenic 
plasma is 



Pr = 



0.01. 



X 



(1) 



where v and x are the diffusion coefficients for momentum 
(by ions only) and thermal energy (by electrons only), re- 
spectively. Pr ~ 0.01 is relatively small, and thus the effects 
of viscosity were expected to be small in comparison with 
those of conduction^ However, the Reynolds number is fairly 
small in the ICM: 
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and thus it i s not s o clea r that the effects of viscosity can 
be neglected. iKunj l|201ll ) (hereafter Kll) recently extended 
the linear dispersion relation for the MTI and HBI to include 
anisotropic viscosity and provides an intuitive, physical expli- 
cation of its effects. The most important is that the growth 
rates of the HBI can be suppressed by viscosity in the limit of 
very rapid conduction (and thus very rapid viscous damping). 
Isotropic viscosity has been ut ilized in a small numb er of pre- 
vious numerical studies; e.g., iRevnolds et~aH (|2005l ) studied 



1 This quantity should not be confused with the often-discussed 
magnetic Prandtl number in which the thermal diffusivity is re- 
placed by the electrical resistivity. The magnetic Prandtl number 
is very large in the ICM. 

2 The effective value of the Prandtl number that enters the per- 
turbed total energy equation in the linear analysis corresponds to 
Pr e g = 0.02. This is because fi = 0.5: only electrons participate in 
conduction, while both electrons and ions contribute to the total 
thermal energy. The net effect is that \ is smaller for the MHD 
fluid by a factor of 2. In our calculations we take Pr = 0.01 and 
At = 0.5. 



the effect of viscosity on the shapes of risin g AGN-blown bub- 
bles. More recently, iDong fc Stond (|2009j) showed that it is 
critical to consider anisotropic viscosity, rather than isotropic 
viscosity, in such calculations because the former is much less 
effective at suppressing the Rayleigh- Taylor instability. 

In this paper, we present fully self-consistent 2D and 3D 
Braginskii-MHD simulations of the ICM, focusing on the evo- 
lution of the MTI and HBI. We introduce our computational 
methods in Sj3] In §4] we present local 2D and 3D simula- 
tions of the HBI and MTI to provide physical insight into 
the role of viscosity. We then use global calculations to study 
the effect of viscosity on the MTI in the outskirts of galaxy 
clusters and the role of the viscous HBI in cluster cores in §5 
& §6, respectively. In the appendix we describe our numeri- 
cal method for anisotropic viscous transport and discuss the 
numerical verification of this algorithm. 



2 METHOD AND MODELS 

We solve the usual equations of magnetohydrodynamics 
(MHD) with the addition of anisotropic thermal conduction 
and anisotropic viscous transport. The MHD equations in 
conservative form are 
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which are the equations of conservation of mass, momentum, 
and energy and the induction equation, respectively. The total 
energy E is given by 



E = e + p 



vv B B 



8tt 



(7) 



where e = p/(j — 1). Throughout this paper, we assume 7 = 
5/3. The anisotropic electron heat flux is given by 

Q = -K Sp bb • VT, (8) 

where ks p is the Spitzer conductivity (Spitzer 1962) and 6 
is a unit vector in the direction of the magnetic field. The 
Spitzer conductivity can be written as «s P = n^sx, where 
X is the actual diffusion coefficient (in units of L 2 T _1 ). The 
viscous stress tensor is given by 



n 



-3pv 



bb-.Vv- 



66- - 
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(9) 



where v is the microphysical momentum diffusion coefficient, 
often termed the kinematic viscosity. Both transport coef- 
ficients are functions of temperature proportional to T 5/l2 , 
where we presume T e = Ti. The microphysics fixes the ra- 
tio of the transport coefficients to be Pr = 0.01 as given by 
Equation (fT]). 

The energy equation also includes a cooling term, C The 
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cooling function we adopt is from iTozzi fc Normanl |200ll ) 
with the functional form 



C = n e n p A(T), 



(10) 



with units of erg cm -3 s _1 . The temperature dependence is a 
fit to cooling dominated by Bremsstrahlung above 1 keV and 
metal lines below 1 keV with 



A(T) = [Ci(k B T) 



- 1 - 7 + C' 2 {k B Tf 5 + C 3 ] 1CT 22 , 



(11) 



where Ci = 8.6 x 10" 3 , C 2 = 5.8 x 10~ 2 , and C 3 = 6.3 x 
10 -2 , for a metallicity of Z = 0.3Zq, with units of [Ci] = 
ergcm 3 s _1 . We use a mean molecular weight of /i ~ 0.62 
which corresponds to a metallicity of approximately 1/3 solar. 

For our simu l ations we use th e Ath ena MHD code 
(|Gardiner fc Stond 120081 ; IStone et al l 120081 ) combined with 
the a nisot ropic conduction methods of IParrish fc Stone] 
(120051 ) and ISharma fc Hammettl |2007l ). The anisotropic vis- 
cosity is implemented in a very similar manner to conduction 
(see the appendix) . The heating, cooling, and anisotropic con- 
duction and viscosity are operator split and sub-cycled with 
respect to the MHD timestep. The cooling simulations are 
implemented with a temperature floor of T = 0.05 keV, be- 
low which UV lines become important, and the cooling curve 
fit is no longer accurate. This temperature floor prevents the 
cooling catastrophe from going to completion. 

This paper will cover a variety of initial conditions from 
local Cartesian boxes to global cluster models. Our initial con- 
ditions will thus be described briefly in subsequent sections 
with appropriate references for more details. In each case, we 
have carried out a least one resolution study to assure numer- 
ical convergence. Any deviations will be noted. 

3 PHYSICS OF THE MTI AND HBI WITH 
VISCOSITY 

We begin with a qualitative description of the physics of the 
HBI and MTI. Despite the mathematical similarities of the 
instabilities, it is useful to discuss each instability separately. 
The MTI is most unstable for horizontal (B _L g) magnetic 
fields with the temperature gradient dT/dz < 0. An upwardly 
displaced fluid element is connected by magnetic field lines to 
a hotter region deeper in the atmosphere. Heat flowing along 
the field expands the fluid element relative to its surround- 
ings, lowers its density and buoyantly destabilizes the per- 
turbation. The upward motion causes the magnetic field to 
be more aligned with the background temperature gradient, 
leading to an instability. 

The HBI, on the other hand, is most unstable for ver- 
tical (B || g) magnetic fields with the temperature gradient 
dT/dz > 0. Imagine a small displacement of fluid elements 
with a wavevector that has a component both parallel and 
perpendicular to gr avity (and B). This configuration, illus- 
trated in Figure 1 of lQuataertl (2008), has regions in which the 
magnetic field lines bunch together and spread apart, leading 
to a converging and diverging heat flux. Rapid heat conduc- 
tion along the perturbed magnetic field lines causes an up- 
wardly displaced fluid element to be heated (by tapping into 
the background heat flux), leading to a buoyant runaway. 

Kll carried out a full linear perturbation analysis on 
equations [3HH] and presents a dispersion relation for the MTI 
and HBI including the effects of anisotropic viscosity. For clar- 
ity when comparing with our results, we reproduce the Kll 



dispersion relation here in a slightly different coordinate sys- 
tem. We define gravity to be in the z-direction, the initial 
magnetic field is in the x-y plane and k\ = fc 2 + fc 2 is the 
component of the wavevector perpendicular to gravity. The 
dispersion relation for the growth rate, a, is given by 



-(6.*) 

IT 2 |ct 2 (a + LO cond ) + aN 2 ^- + CJcondCJbuoy/cJ 

a \a 2 (a + cj cond ) + (aN 2 + ^cond^b uoy ) j 



(12) 
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where a 2 = a 2 -(k ■ va) 2 , and va = B/(47rp) 1/2 is the Alfven 
velocity. The Brunt- Vaisala, oscillation frequency is given by 

, r2 1 dPdlnS 



dz 



(13) 



7P dz 

where S = Pp 1 , and corresponds to buoyant oscillations in 
an unmagnetized plasma (g-modes). The characteristic fre- 
quency at which conduction and viscosity act are given by 

2^ 
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where \ ls the thermal diffusivity and va is the ion collision 
frequency. The buoyancy frequency is given by 
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dz 
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which is roughly the fastest growth rate of the MTI or HBI. 
Finally, the dimensionless geometric factor is given by 



K, 



2bi 



+ ■ 



A- 2 



(17) 



where b x — B x /B is the dimensionless magnetic field 
strength. 

We can greatly simplify this expression and improve our 
physical intuition by assuming that the magnetic field is weak 
(or, equivalently, that the perturbation wavelength in the 
cluster core A>«a/t~2 kpc, where r is a typical timescale 
for the instabilities to grow). We also assume the perturba- 
tion lies in the plane containing B and g, i.e. k y — 0. In this 
limit, the dispersion relation simplifies to 



a (a + ^cond) 



b k 



(18) 



+ O'A^-p- + Wcond^buoy'C = 0. 



Without viscosity, there is a fast conduction limit in 
which the growth rate asymptotes to a constant a — 
Wbuoyfex/fc 2 ; with viscosity, however, no such limit exists. 
In Figure [Tj we plot the theoretical curves for the growth 
rates of both the HBI and MTI as a function of the ratio 
^cond/i^buoy for several Prandtl numbers. These curves as- 
sume a wave vector k at 45 degrees from the magnetic field 
B; these represent typical, not maximum, growth rates. The 
ratio u con d/ti)buoy oc k 2 so that larger values of this ratio cor- 
respond physically to smaller scales at fixed conductivity. 

Figure [Tj shows that the inviscid solution reaches an 
asymptotic growth rate in the fast conduction limit; how- 
ever, with viscosity the growth rate slowly decreases as one 
moves to the infinite conduction limit. For the HBI, the 
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Figure 1. The lines show the theoretical linear growth rates of the MTI (right) and HBI (left) for a variety of conduction frequencies and 
several Prandtl numbers for k at 45° relative to the magnetic field. This geometry is just an example and does not show the fastest growing 
mode for either the MTI or HBI. In particular, for the MTI there are modes with growth rates comparable to the inviscid (Pr = 0) case 
even for w con( j 3> Ubuoy (Kll); this is not true for the HBI. Real astrophysical plasmas have Pr = 0.01. Pr = 0.06 is shown only as a test 
of our numerical methods. Low w con< j corresponds to larger scales at fixed conductivity. Measured growth rates from 2D simulations are 
plotted as star symbols and agree well with the analytic results. 



fastest growth occurs on scales satisfying o; v i sc ~ uibuoy, i.e., 
kJcond ~ 6ojbuoy, and the maximum growth rate for Pr = 0.01 
is about 70% of the asymptotic inviscid value. This is in- 
dependent of k± unless k± <C fen . For higher viscosities, 
the maximum growth rate scales as ~ t^buoy /Wise • The in- 
fluence of viscosity on the MTI is much more modest. In 
particular, the fastest growing mode is achievable even for 
tj v i sc , aj C ond S> kJbuo y ; this is not seen in Figure [T] because of 
the particular perturbation chosen. One can, however, always 
find short wavelength MTI modes that grow at ~ £Jb UO y 

The microphysical plasma physics fixes the ratio of the 
diffusion coefficients, the Prandtl number (eqn. [I]), to be 
Pr = 0.01. This value corresponds to o; v i sc /aj con d = 1/6. 
Because there are MTI modes that grow rapidly even for 
^cond 2> ^buoy , we would a priori not expect viscosity to have 
a significant effect on the evolution of the MTI. For the HBI, 
the viscous and inviscid growth rates are similar only when 
Woond/^buoy ^ 10; these larger scales in a physical system 
are only minimally affected by viscosity. We thus expect the 
HBI to be modified by viscosity only if cj con d 2> 10 ^buoy on 
the largest scales of the system. 

In Figure[2]we show the ratio of Wcond/wbuoy f° r observed 
clusters using k — 2n/r as the conduction length scale and 
not including any geometric factors, i.e. k ■ b — 1. The left 
panel focuses on the core of the cluster where the HBI op- 
erates while the right panel is at larger radii where the MTI 
is present. The values of w oon d/Wuoy m Figure [2] are charac- 
teristic of the largest scales in the system. Smaller scales are 
more viscous/conducting, but are also where magnetic ten- 
sion is most likely to suppress the MTI/HBI. Figure [2] also 
shows the ratio o; con d/ajbuoy (defined in the same way) for our 
model clusters used in the global simulations in §5 & 6. The 
models are generally representative of real clusters. 

The data in Figur e [5] comes from the ACCEPT sample 
(|Cavagnolo et al . 2009); the analysis of the data follows the 



methodology described in iMcCourt etHI l|2011bl ) with the 
clusters listed in Table 2 of that work. These examples span 
relatively massive NCC clusters all the way down to groups. 
Cool-core (CC) clusters have u con d/uJbuoy ~ 5-30 in the bulk 
of the core between a few-200 kpc. By comparing the real 
CC clusters to Figure [1] we see that viscosity only modestly 
reduces the growth rates of the HBI. Non-cool-core (NCC) 
clusters are more conducting/viscous as a result of the hotter 
and lower density cores. The effects of viscosity on the HBI 
are theoretically more likely to be significant in NCC clusters. 

The star symbols in in Figure [T] show good agreement 
between the measured growth rates in our simulations and the 
analytic results; this represents a strong test of our algorithms 
for anisotropic conduction and viscosity. We discuss this in 
more detail in the appendix. 



4 LOCAL SIMULATIONS 
4.1 Initial Conditions 

It is illustrative to start with the simplest possible experi- 
ments, namely local two- and three-dimensional boxes with 
system sizes L < the scale-height H. For this section we work 
in units with fcg = m p = 1, go — — 1, and with a hydrogenic 
plasma t hat has \i = 1/2 . Our i nitial conditions are fully de- 
tailed in IMcCourt etaLl (|2011al ). For the HBI, we start with 
a simple hydrostatic equilibrium that is linearly unstable: 

T{z) = T (l + z/H), (19) 
p(z) = Pa {l + z/H)-\ (20) 

where we set To = po — 1 and choose a scale height of H = 2. 
Our boxes are of size L = 0.2 in each dimension with a resolu- 
tion of (96) 2 or (96) 3 . Results at this resolution are very well 
converged by all metrics. The vertical boundary conditions 
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Figure 2. Left: The measured ratio of w cont i/kJbuoy f° r the cores of clusters in the ACCEPT database, taking an average value of 
|dlnT/dlnr| = 1/3 in t^buoy and k = 2n/r in oj con( j. The latter choice corresponds to the conduction time across the local radius r and 
represents the smallest value of w con< j/o;i, UO y at a given radius. Cool-core clusters have ^cond/^buoy ~ 10 m the bulk of the cool core. 
Non-cool-core clusters have somewhat more rapid conduction and are thus more viscous as well. We have also overplotted the values of 
w cond/ w buoy m our global cluster models used for HBI (the CC and NCC HBI models). Right: The measured values of t^cond/^buoy m the 
outskirts of clusters in the ACCEPT sample. Our model cluster for the global MTI simulation is overplotted as well. Our fiducial model 
clusters are reasonably similar to observed clusters. 



fix the temperature to its initial value, and the horizontal 
boundary conditions are periodic. The initial magnetic field 
is weak and vertical, -Bo = 10 _6 z. 

Fo r the local MTI simulat ions we utilize the set-up from 
§3.3 of iMcCourt et all (|2011ah given by 



TO) = To exp 



S^buoy 

90 



(l-e^ s ) 



g(z) = g e~ 



z/S 



with To = go = 1 and S = 3. We select lu'^ u 



(21) 

(22) 
: 1/2 and 

solve numerically for p(z) to ensure hydrostatic equilibrium. 
These boxes are of size H/2 with an initially weak horizontal 
magnetic field, So = 10 _6 z. Convergence is a bit more subtle 
in these simulations, b ut 128 3 is reasonably well-converged 
(|McCourt et alj|2011al ). 

For the HBI, the atmosphere satisfies the Schwarzschild 
criterion (dS/dz > 0) and would be buoyantly stable in the 
absence of anisotropic conduction. In these calculations we fix 
the difiusivities so that the conduction frequency across the 
box is a chosen constant relative to the buoyancy frequency, 
and show results for different values of w coll( j/cJbuoy for both 
the inviscid case and the physical case with Pr = 0.01. 



4.2 Nonlinear Saturation 

4.2.1 HBI 

We begin by examining the nonlinear saturation of the HBI 
with 2D single mode simulations. The magnetic field lines 
from these simulations are visualized in Figure Each verti- 
cal column compares inviscid and viscous simulations at the 



time t = 30 £hbi • This is well into the non- linear regime, in 
which the statistical properties of the plasma do not evolve 
strongly with time. We are able to simulate a variety of phys- 
ical scales in these local simulations by changing the constant 
conductivity. We then label the simulation with a ratio of 
i^cond /^buoy calculated using the box size, A = L, for the 
conduction length in cj C ond- For w C ond/i^buoy = 2, the left- 
most column shows that viscosity does not strongly change 
the saturated state. The inviscid cases shown in the top row 
all show qualitatively no difference in the nonlinear satu- 
rated state as a function of cj con d/o;buoy; however, the ef- 
fect of anisotropic viscosity becomes clear for the viscous case 
with w C ond A-'buoy = 60 (bottom right). This case shows promi- 
nent, highly-bunched vertical magnetic field structures that 
are never seen in the inviscid case. This bunching occurs be- 
cause the field lines cannot slip past each other in 2D, a point 
we will return to shortly. 

In order to better understand the limitations of 2D sim- 
ulations, it is useful to compare the 2D and 3D evolution for 
identical parameters. Figure [4] shows the exact same simula- 
tions performed in 2D and 3D for three different conductivi- 
ties. For the HBI, the results for simulations with different box 
sizes are essentially the same pro vided that the same va lue of 
^cond (k = 27r/L)/w b uoy is used l|McCourt et al. l l2011al ). We 
demonstrate this result explicitly in Figure [S] which compares 
the magnetic geometry evolution for the viscous HBI in two 
different box sizes. To keep the ratio of Wcond/^buoy = 50 
fixed, the simulation with L/H = 0.3 has a conductivity and 
viscosity that are 9 times larger than the simulation with 
L/H = 0.1. The evolution, especially at linear and early non- 
linear times, is nearly indistinguishable. At very late times, 
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Figure 3. Snapshots of the magnetic field lines from 2D HBI simulations. The top row displays inviscid simulations, while the bottom row 
includes anisotropic viscosity. The conductivity increases from left to right; the snapshots at a given conductivity are shown at the same 
absolute time for both the inviscid and viscous simulation (30 HBI growth times for the inviscid model). At low ^cond/o^oy, viscosity has 
little effect on the nonlinear evolution of the HBI; however, for ^JcondA^buoy ~ 60, anisotropic viscosity causes the magnetic field to bunch 
into pronounced vertical structures. This effect is much weaker in 3D simulations relative to the 2D results shown here (see Fig. [4}. 



the larger box has experienced slightly more magnetic field 
evolution as a result of the different modes in the domain 
relative to the scale height. We conclude from Figure [5] that 
local simulations of the HBI with L < H are independent 
of domain size; the results depend instead on the value of 
^cond/^buoy, the fundamental dimensionless parameter that 
quantifies the influence of conduction and viscosity. Thus al- 
though the boxes in Figure [3] are only ~ 0.1 H in size, the 
local simulations with w con d/wbuoy = 18 & 50 are quite in- 
dicative of how the HBI would evolve on large scales in the 
ICM (see Fig.[5|. The large scales are the most important be- 
cause both tension and viscosity suppress the growth of the 
HBI on smaller scales. We justify this interpretation of the 
local HBI simulations using global cluster models in §5] 

Figure U shows the mean magnetic field direction, (|6 Z |), 
which is related to the effective vertical thermal conductivity 
as fs P ~ {b z }- In 2D the ability of the HBI to reorient that 
magnetic field is retarded by viscosity relative to the invis- 
cid case, especially for aj C ond/^buoy = 50 (right). However, in 
3D the effect of the viscosity is much less pronounced. The 
biggest effect of viscosity is that for tJ con d/u;buoy = 50, the 
initial growth of the HBI is slower; however, the instability 
still produces a significant rearrangement of the magnetic field 
by t ~ 50 tbuoy • We note that our initial velocity perturba- 
tions have white noise amplitudes with a mean magnitude of 
5 x 10~ 4 c s . For a larger and more realistic perturbation, the 
linear phase would complete more quickly. 



What is the cause of the striking difference in the 2D 
and 3D evolution of the HBI with viscosity? Recall in Fig- 
ure [3] that the magnetic field lines strongly bunched up in 
2D and could not slip past each other. This type of motion 
is known as the interchange mode and can be thought of as 
two magnetic flux tubes slipping past each other. Interchange 
motions have »1B and therefore are not damped at all by 
anisotropic viscosity, although they would be damped by an 
isotropic viscosity. In 3D the magnetic field lines are able to 
slip past each other, and the HBI is able to continue to re- 
orient the magnetic field. This effect is seen in a different 
guise i n studies of the magnet ic Ray leigh- Taylor (RT) insta- 
bility jStone fc Gardiner|[2007h . In 2D, magnetic tension sup- 
presses the RT instability below a critical wavelength; how- 
ever, in 3D magnetic tension does nothing to suppress the 
interchange motions of flux tubes. Thus, the instability grows 
despite the apparent wavelength cut-off in the dispersion re- 
lation. This uniquely 3D effect of anisotropic viscosity makes 
2D studies of nonlinear saturation irrelevant, and greatly re- 
duces the effect of viscosity relative to what one might naively 
expect from the dispersion relation. It is particularly striking 
for aj C ond/f-'buoy = 50 (right panel in Fig. [4} the non-linear 
evolution in 3D is quite similar with and without viscosity, 
although the viscous evolution is delayed by the considerably 
longer linear growth time. 
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local simulations (L = 0.1H). We quantify the field evolution using the vertical component of the magnetic field. The differences between 
the viscous and inviscid simulations are much smaller in 3D (blue lines) than in 2D (black lines). We attribute this to the presence of 
interchange modes in 3D. Even for ^cond/^buoy ~ 50 the nonlinear evolution of the HBI in 3D is similar to that in inviscid simulations, 
with the primary difference being that the nonlinear satuation is delayed because of the reduced linear growth rate. 



Figure 4. The nonlinear evolution of the magnetic field due to the HBI for £J con d/tUb u 



4.2.2 MTI 

Next, we also examine the nonlinear evolution of the MTI in 
3D Cartesian simulations. The dispersion relation shows that 
the fastest growing modes of the MTI are less affected by vis- 
cosity than they are for the HBI. Figure [6] shows the results 
of our Cartesian simulations of the MTI with and without 
viscosity and reveals several interesting properties. First, the 
saturated kinetic energies are largely independent of both the 
ratio of w C ond/ w buoy an d whether anisotropic viscosity is in- 
cluded. Note that the simulation with w con d /c^buoy = 50 is 
consistent with the range of clusters in the ACCEPT sam- 
ple that we have examined (see the right panel of Figure 
The fact that the kinetic energies are so similar even in the 
high-viscosity simulation, is likely the result of interchange- 
like turbulence. 

Second, we find the very interesting result that the 
growth of the magnetic energy (middle panel of Figure 
is suppressed as the viscosity is increased. This result can be 
understood in relatively easy terms: the increase in the mag- 
netic field strength in this type of turbulence is proportional 
to the increase in the length of the magnetic field line, as it 
is stretched and tangled. These parallel stretching motions 
are precisely those motions that are damped by anisotropic 
viscosity, thus suppressing the growth of the magnetic field. 

Finally, we consider the reorientation of the magnetic ge- 
ometry by the MTI (bottom panel of Figure [6]). In the invis- 
cid case, an initially horizontal magnetic field is re-oriented 
to be largely isotropic, (|6 Z |) = 1/2 in 3D. For modest vis- 
cosities ^ con d/a;buoy = 1, there is little change in the sat- 
urated geometry; however, for higher conductivity and vis- 
cosities (e.g., cj con( j /c^buoy = 50) the magnetic field becomes 



slightly more vertical than isotropic (about a 10% change 
in (\b z \)). We can gain insight into this behavior by consid- 
ering the evolution of simulations with an initially vertical 
field, a case tha t is linearly stable, bu t non-linearly unstable 
(see Figure 6 of Mc Court et aljfefjllal ). An initial horizontal 
perturbation is purely Alfvenic in nature, and thus is unaf- 
fected by Braginskii viscosity. A simulation with an initially 
vertical field and ^cond/^buoy = 1 (blue dotted line) evolves 
towards an isotropic magnetic field configuration. However, 
in doing so, the field develops a component parallel to the 
velocity, which is susceptible to viscous damping. Thus, for 
an initially vertical magnetic field, the higher viscosity sim- 
ulation (^cond/wbuoy = 50, blue dashed line) evolves much 
more slowly than its less viscous counterpart. It is precisely 
this partial stabilization of modes that initially have the field 
aligned with gravity that results in the radial bias seen in the 
high-conductivity, horizontal simulation. 



5 GLOBAL MODELS OF THE MTI IN 
CLUSTER OUTSKIRTS 

In this section we discuss the effect of Braginskii viscosity 
on the evolution of the MTI in the outskirts of galaxy clus- 
ters using global cluster models. Beyond approximately the 
scale radius of a cluster, the temperature profile of the ICM 
almost always declines with radius, thus making it unstable 
to the MTI. Recent work has shown that the MTI is able 
to drive large turbulent velocities which provide non-therma l 
pressure support in hot, massive clusters (|Parrish et alj201lh . 
It is critical to understand this non-thermal pressure support 
and its implications for measuring cluster masses for use in 
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Figure 5. The evolution of the magnetic field orientation for the 
viscous HBI (Pr = 0.01) for two different domain sizes relative to 
the temperature scale height (L/H). To keep the ratio of w con( j(A: = 
2tt/L) /cjbuoy = 50 fixed, the simulation with L/H = 0.3 (dotted 
line) has a conductivity and viscosity that are 9 times larger than 
the simulation with L/H = 0.1 (solid line). For local simulations 
(L < H) the linear and nonlinear evolution are nearly independent 

of box size for fixed tJ C ond A^buoy ■ 



cosmological parameter estimation through either X-ray or 
Sunyaev-Zelovich (SZ) methods. 

Our initial condition for this section is a spherically- 
symmetric, hot, massive cluster that resembles Abell 1576 
with a mass of 1.6 x 10 15 Mq. We use a softened NFW pro- 
file with a scale radius of r s = 600 kpc and a softening ra- 
dius of 70 kpc. We initialize an atmosphere in hydrostatic 
equilibrium using the entropy power law in the ACCEPT 
database for Abell 1576: a central entropy Ko = 186 keVcm 2 , 
K\ = 98 keVcm 2 , and power-law exponent, a — 1.38 
(jCavagnolo et al1l2009l ) . If we assume that our fiducial cluster 
is located at z — 0.1, then for the appropriate WMAP5 cos- 
mological parameters rsoo = 1-09 Mpc, and the virial radius 
is V2oo ~ 1.6 Mpc, where ta corresponds to an overdensity 
of A times the critical density. We do not include cooling, as 
we focus on the portion of the cluster that is w ell outside the 
coolin g radius. For full simulation details, see IParrish et alJ 
(|20T lh . By examining Figure |5J we can see that the ratio of 
Wcond/wbuoy m our m °del is consistent with the range ob- 
served in real clusters in the ACCEPT database at the radii 
of interest (r > 200 kpc). 

The simulations are carried out on a (196) 3 Cartesian 
grid in a computational domain that extends from the center 
of the cluster out to ±1300 kpc. Within this Cartesian do- 
main, we define a spherical subvolume with a radius of 1225 
kpc from which we extract cluster properties. In this volume 
we initialize tangled magnetic fields with (\B\) = 10~ 8 G 
(plasma ft ~ 10 4 -10 6 ) and a Kolmogorov power spectrum. In 
order to simulate clusters with a negative radial temperature 




buoy 



Figure 6. The nonlinear evolution of the MTI for w con d /u>buoy = 1 
and 50 in 3D (L = H/2) with and without viscosity. We omit the 
u cond/ u buoy = 50 inviscid simulation, as its evolution is almost 
exactly the same as the w con( j /i^buoy = 1 inviscid simulation. The 
kinetic energy (top panel) is largely independent of viscosity. The 
magnetic energy generation (middle panel) is suppressed for higher 
viscosities. The magnetic geometry (bottom panel) for initially hor- 
izontal magnetic fields (black lines) becomes relatively isotropic for 
all the simulations, although slightly more vertical with higher vis- 
cosity. Statistical isotropy corresponds to (|6 Z |) = 0.5 (red, long- 
dashed line). Simulations with initially ve rtical fields (blue lines ) 
are linearly stable but nonlinearly unstable ( McCourt et al . 2011a); 
viscosity slows the reorientation of magnetic field in this case. 



gradient, we fix the temperature at the peak of the cluster 
temperature profile (approximately 10 keV at 200 kpc) and 
at the maximum radius (approximately 6.5 keV at 1225 kpc) 
of our model cluster to the initially-computed temperature 
values. Thus, we are imposing Dirichlet boundary conditions 
on the temperature profile. This fixed temperature gradient 
then drives the continued evolution of the MTI. 

Figure[7]shows the non-linear evolution of the MTI in the 
global galaxy cluster context. The evolution of the kinetic 
energy with time is quite similar for both the viscous and 
inviscid case just as in the local simulations (Figure [B); how- 
ever, the amplification of the magnetic field is again somewhat 
suppressed in the viscous case. The magnetic field stretching 
necessary to amplify the field is exactly the motion that is 
damped by anisotropic viscosity. This damping of magnetic 
field amplification makes it more difficult to use a turbulent 
dynamo, regardless of the source of the turbulence, to amplify 
the magnetic field from primordial values to those observed 
today in galaxy clusters. The magnetic geometry evolution is 
very similar with and without viscosity, as we initially began 
with a geometrically isotropic magnetic field. A very slight 
radial bias exists in the steady state in both cases. 

Figure [S] shows the saturated and azimuthally-averaged 
rms Mach number profile for our model cluster at a time of 
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Figure 7. The non-linear evolution of volume-averaged properties 
of the MTI in clusters in our global model with (dotted lines) and 
without (solid line) viscosity. The evolution of the volume-averaged 
kinetic energy (top panel) and the the magnetic field geometry 
(bottom panel) are largely unchanged by the viscosity. The mag- 
netic geometry is measured with respect to the radial direction: 
9 = 0° corresponds to radial, and 6 = 60° (red, long-dashed line) 
corresponds to a random field. The amplification of the volume- 
averaged magnetic energy (middle panel) is modestly suppressed 
by anisotropic viscosity, similar to what was observed in our local 
simulations. 



8.3 Gyr. The turbulence profiles with and without viscosity 
are almost identical beyond 400 kpc. More quantitatively, at 
rsoo the rms Mach numbers differ by only 5%. Near the very 
center of the cluster, the rms velocity is actually larger for the 
viscous simulation. Although we have no definitive explana- 
tion, this is perhaps due to the turbulent eddies becoming less 
aligned with the magnetic geometry. Since our cluster is on 
the massive and hot end of the cluster mass function, this clus- 
ter has a particularly large value of ^cond/wbuoy, which yields 
the greatest effect of viscosity on the linear MTI (see Fig- 
ure [T]). Since viscosity has little effect on the MTI-generated 
turbulence in this hot cluster, we expect viscosity to make 
almost no difference for the MTI in lower mass (and cooler) 
clusters. 



6 GLOBAL MODELS OF THE HBI IN 
CLUSTER CORES 

We now move inwards to consider the effects of viscosity in 
global models of the cores of galaxy clusters. Within 200 kpc, 
the cooling times are as short as 100 Myr at ~ 10 kpc. These 
short cooling times, along with evidence that large mass fluxes 
of gas are not cooling from a ho t to cold phase constitut e 
the modern cooling flow problem (|Peterson fc Fab ian 2006). 
Namely, what is keeping cool core cluster cools from cooling 
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Figure 8. The saturated and azimuthally-averaged rms Mach 
number profile due to the MTI in our global model cluster (model 
MTI in Fig. [2j. The inviscid (solid line) and viscous (dotted line) 
simulations differ only by a small amount. 



further? There has been speculation that thermal conduc- 
tion from large radii could o ffset the cooling luminosity (e.g., 
iNaravan fc Medvedevll2~00ll ): however, the HBI presents a se- 
rious problem to this scenario. Alternatively, there are many 
observations of AGN feedback that present a more plausible 
heating mechanism. 

The first cluster model we consider is m odeled on obser- 
vations of Abell 2199 (| Johnstone et al.ll2002T i. We initialize a 
cluster in a static NFW halo with a mass of 3.8 x 10 14 M Q 
and a scale radius of 390 kpc. In this potential we calcu- 
late a spherically-symmetric atmosphere in hydrostatic equi- 
librium and thermal equilibrium with conduction (at 1/3 of 
Spitzer) exactly balancing cooling. The model cluster has a 
central temperature and electron density of ~ 2.0 keV and 
~ 0.021 cm~ 3 , respectively, and a temperature and density of 
5 keV and 1.67 x 10~ 3 cm" 3 at 200 kpc. This corresponds to 
a central cooling time of 1.7 Gyr. In this region we initialize 
a tangled magnetic field with a Kolmogorov power spectrum 
and an amplitude of (\B\) ~ 10~ 8 . The magnetic field is tan- 
gled on scales from 50 kpc down to 30 kpc. The simulations 
are computed on a Cartesian grid in a domain extending from 
the cluster center to 240 kpc with (128) 3 gridpoints. A resolu- 
tion study at (256) 3 showed that the behavior of the magnetic 
field geometry was almost indistinguishable from the lower 
resolution results we focus on here (Fig. [9]). The magnetic 
field geometry and temperature profiles as a function of ra- 
dius (not plotted) are nearly identical; as a result, we consider 
these simulations very well-converged. In these simulations, 
the temperature is fixed to the initial temperature at a ra- 
dius of 200 kpc everywhere, but there is no cen tral boundary 
condit ion. For more details of the set-up see iParrish et al.l 
(2009). Our CC atmosphere model has u con d /^tmoy ~ 10, 
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Figure 9. Time evolution of the volume-averaged angle of the mag- 
netic field with respect to the radial direction in our global cool-core 
cluster model (model CC in Fig.rjl. 9 = 0° is radial. 9 = 60° cor- 
responds to a random, isotropic magnetic field. The HBI reorients 
the magnetic field, reducing the effective conductivity, which pre- 
cipitates the cooling catastrophe at 2.4 Gyr. Very little difference 
is seen between the inviscid (solid line) and viscous (dotted line) 
evolution. We demonstrate that these results are well-converged by 
plotting the nearly identical evolution of a viscous simulation with 
double the resolution (red dashed line). 
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Figure 10. Azimuthally-averaged temperature profiles for iden- 
tical cluster cores with different imposed rms turbulent velocities 
(see legend) with a fixed driving scale of L = 40 kpc. A very strong 
bimodality is seen in the stability properties of the thermal profiles. 
Braginskii viscosity is included. Driving at 71 .7 km s — 1 produces a 
stable temperature profile; whereas a tiny change in the driving to 
70.5 km s — 1 leads to the cooling catastrophe. The two thermally 
stable runs report the measured time in parentheses; while the two 
lower velocity unstable runs (marked with a *) report the time of 
cooling catastrophe in parentheses. 



which is reasonably consistent with the cool-core ACCEPT 
cluster plotted in Figure [2] Relative to observed clusters, our 
model is somewhat high in ^cond/^buoy at small radii and on 
the low-end at larger radii. This discrepancy is due to the 
fact that by enforcing thermal equilibrium, we demand that 
conduction match cooling everywhere without any AGN feed- 
back. This is, perhaps, a hint that other heating mechanisms, 
such as AGN feedback, are necessary for real CC clusters. 

We first evolve our fiducial cluster model with anisotropic 
conduction and cooling only. We diagnose the effect of the 
HBI by calculating the mean magnetic field angle from radial 
(8) — cos _1 (|b ■ r\). This quantity is relevant to the thermal 
evolution as the effective radial thermal conductivity, often 
called the Spitzer fraction, is given by 

/ S p=Qr-/Qr~COS 2 <|S-f|), (23) 

where Q r = — KdT/dr, which is the radial heat flux if the 
conduction were purely isotropic at the Spitzer value. Figure 
[9] shows that the HBI acts to reorient the magnetic field to be 
more azimuthal, thus reducing the effective radial conductiv- 
ity. In this calculation the cooling catastrophe (defined as the 
central temperature reaching our imposed floor of 0.05 keV) 
occurs at ~ 2.4 Gyr. 

When the effects of Braginskii viscosity are added, the 
evolution changes very little as shown by the dotted line in 
Figure As we saw previously in the local calculations, in 
3D interchange-type motions are able to proceed unimpeded 



by any viscosity at all. As the HBI drives rather small turbu- 
lent velocities (12 km s _1 at 2 Gyr), the viscous dissipation 
does next to nothing to impede the cooling catastrophe which 
occurs at almost the exact same time as in the inviscid simula- 
tion. In the Braginskii (collisional, anisotropic) limit, heating 
by viscous dissipation of turbulence is thermally unstable as 
V • (n • v) oc v oc T 5 / 2 i go even more vigorous turbulence 
cannot stably balance cooling. 

Turbulence driven by galaxy wakes, AGN feedback, 
structure formation, or any other source can counteract 
the ability of the H BI to reorient magnetic field lines 
(|McCourt et al.ll2011al ). Quantitatively, turbulence on a scale 
L with velocity Sv is able to suppress the HBI when 

L ^ + J dlnTV 1/2 

Wd y (L) ~ — < CteBi k i I g-^- J , (24) 

where tuBi is the HBI growth time, and f is a dimension- 
less constant that is deter mined by simulations (I Sharma et alj 
120091 ; IParrish et alJlioibT ). Since t cddy oc L 2/3 , this inequality 
is most difficult to satisfy at the outer scale, which domi- 
nates the turb ulent energy. Pos sible sources of turbulence are 
galaxy wakes |Kim et al.ll2005h or AGN feedback itself. 

We simulate the interaction of turbulence, the HBI, and 
cooling in the fiducial cluster core model described previously 
by adding a random velocity forcing. We drive the velocity 
fields in Fourier space with a flat spectrum on a scale L = 
40 ± 10 kpc such that the phases are random and the forcing 
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is incompressible. We have simulated forcing that results in 
turbulent velocities from 30-120 km s _1 , which corresponds 
to a pressure support of less than a few percent in the cluster 
cores. The turbulent energy input from this forcing. e tur b — 
p(5v) 3 /L, is negligible compared to the cooling rate, and thus 
the turbulence is not energetically important. 

Figure mJl shows the late-time, azimuthally-averaged tem- 
perature profiles for the same cluster model with differ- 
ent injected turbulent velocities. Viscosity is present at the 
Spitzer value for these simulations. At low turbulent veloc- 
ities, 5v < 70.5 km s _1 , the HBI acts to reorient the mag- 
netic field lines, reducing the effective thermal conductivity 
and driving the cluster core towards a cooling catastrophe. 
The time of the cooling catastrophe is reported in parenthe- 
ses in the legend. Turbulence does indeed delay the cooling 
catastrophe somewhat, but below a critical value is incapable 
of preventing it. For example, a turbulent forcing velocity of 
70.5 km s _1 is able to delay the cooling catastrophe from 2.4 
Gyr with no turbulence to 12 Gyr with turbulence. 

On the other hand, for larger turbulent velocities, the 
turbulence is able to keep the magnetic field geometry close 
to isotropic. In this case, the central temperature approaches 
a new thermally stable state at a higher temperature. In fact, 
as 8v increases, the central temperature increases due to both 
conduction and the larger radially-inward advective heat flux 
driven by the turbulence. In the final state of the cluster, 
the central entropy is significantly larger than in the initial 
condition. For example, a driving velocity of Sv ~ 86.0 km s~ 
takes the initial central entropy from Ko ~ 26.5 keVcm 2 to 
90 keVcm 2 . 

We find a strikingly strong bimodality in the behavior 
of the cluster temperature profiles and central entropies. Fig- 
ure \TU\ shows that changing the driving velocity by a mere 
1.2 km s _1 causes a cluster to go from heating approximately 
balancing cooling (71.7 km s _1 ) to a core experiencing a cool- 
ing castatrophe (70.5 km s _1 ). The former simulation is al- 
most exactly on the border of the stability boundary and 
departs only marginally from the initial temperature profile. 
This particular case where conduction plus turbulence is bal- 
ancing cooling to maintain a cool core is an example of fine 
tuning. For turbulent velocities 5% above or below this value, 
the core evolves to a clear non-cool-core (NCC) profile or 
experiences a cooling catastrophe. Note, that we have not 
included any source of AGN feedback that could avert the 
cooling catastrophe in the latter situation. 

What does Braginskii viscosity do? Including the viscos- 
ity makes little difference to the fundamental bimodality ob- 
served. In fact, it may even sharpen the transition since the 
viscosity scales as T 5 / 2 . A cluster whose core is heating from 
the initial condition becomes more viscous and thus a given 
turbulent velocity is more readily able to reorient the mag- 
netic field relative to the HBI as the temperature increases. 

We checked the robustness of our conclusions by varying 
our initial fiducial physical a variety of ways: 

• A hotter, more massive CC cluster. We initialize a model 
with a central temperature of 2.5 keV and entropy of Ko = 
27 keV cm 2 which increases to a temperature of 8 keV at 300 
kpc. In the absence of turbulence this cluster proceeds to a 
cooling catastrophe in ~ 2.1 Gyr, much as our fiducial model. 
With turbulent driving at an rms velocity of 57 km s _1 , the 
cluster transitions to a NCC cluster with a central entropy 



of Ko = 71 keV cm 2 . The viscous and inviscid cases are very 
similar as in our fiducial calculation. 

• A NCC cluster of equivalent mass. We initialize our fidu- 
cial model with temperature ranging from 4 keV to 5 keV 
over 200 kpc and a central entropy of 98 keV cm 2 . From linear 
theory, one expects viscosity to reduce the HBI growth rates 
somewhat due to the higher ratio of w C ond/^buoy relative to 
our fiducial model (see Figure [TJ ; however, the evolution of 
the temperature profile is almost indistinguishable with and 
without viscosity. Without any imposed turbulence this clus- 
ter reaches a cooling catastrophe around 6.2 Gyr which is 
somewhat faster than the initially predicted cooling time of 
7.2 Gyr0 

• A hotter, more massive NCC cluster We initialize a hot- 
ter cluster with a mass of 5.2 x 10 M§ and a temperature 
that ranges from 6 keV at the center to 8 keV at 300 kpc. 
This model has a scale radius of 650 kpc and a central en- 
tropy of Ko = 129 keVcm 2 . The profile of oj con d/cJbuoy is 
shown in Figure [2] with the NCC (HBI) label. The value of 
^cond/^buoy is ~ 100 at 50 kpc and ~ 40 at 100 kpc, consis- 
tent with the ACCEPT NCC measurements. Figure [TT] shows 
the evolution of the radial profile of the magnetic field angle 
from its initial statistical isotropy of 60°. In the majority of 
the volume, the magnetic field is reoriented to be substan- 
tially more tangential (8 > 75°) with little difference between 
the viscous and inviscid cases. The small radial bias around 
20 kpc is caused by inhomogeneous radial infall as the cluster 
begins to suffer a cooling catastrophe. Based on the profiles 
in Figure viscosity is even more important at these small 
radii in the initial cluster model; additional suppression of 
the HBI may be present but is overwhelmed by the effects 
of cooling in our global calculations. At larger radii, viscous 
effects are smaller, and the HBI has already reduced conduc- 
tion from the outskirts. Overall this hot, NCC cluster model 
shows that the HBI can operate even when viscous effects are 
important. 

• Reduced magnetic field correlation length. We simulate 
our fiducial cluster with the magnetic field tangled on scales 
ranging from 20 kpc to 10 kpc. The shorter correlation length 
can reduce the effective conduction length scale, increasing 
^cond, and reducing the HBI growth rate. The thermal evo- 
lution, however, is nearly identical to our fiducial model both 
with and without viscosity. 

• Higher initial magnetic field. We simulate our fidu- 
cial cluster with the magnetic field increased to (\B\) ~ 
10~ 6 which is similar to observed magnetic field strengths 
(|Carilli fc Tavlorll2002D . Magnetic tension is able to suppress 
short wavelength modes in a way very similar to viscosity. 
With the higher field, the HBI proceeds slightly more slowly 
and the fiducial cluster reaches a cooling catastrophe approx- 
imately 200 Myr later. We find no difference between the 
viscous and inviscid simulations. 



7 DISCUSSION AND CONCLUSIONS 

The dilute plasma in the ICM of galaxy clusters has a mean 
free path along magnetic field lines that is many orders of 
magnitude larger than the electron/proton gyroradii. As a 

3 This cooling time is an overestimate as the Bremsstrahlung cool- 
ing rate increases oc T — 3 / 2 as material cools at constant pressure. 
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Figure 11. Radial profiles of the azimuthally-averaged angle of the 
magnetic field with respect to the radial direction in our global non- 
cool-core cluster model (model NCC in Fig. [2} at two times. 6> = 0° 
is radial. 9 = 60° corresponds to a random, isotropic magnetic 
field. The HBI reorients the magnetic field, reducing the effective 
conductivity, which precipitates the cooling catastrophe at 2.4 Gyr. 
Very little difference is seen between the inviscid (black curves) and 
viscous (red curves) evolution. The radial bias seen around 20 kpc 
is due to inhomogcncous radial infall. 



result, thermal conduction by electrons and momentum trans- 
port by ions are strongly anisotropic with respect to the mag- 
netic field. In this paper we have carried out the first simula- 
tions of buoyancy instabilities in the ICM (the MTI and HBI) 
that include anisotropic conduction and viscosity simultane- 
ously and self-consistently. This extends the analytic work of 
Kll, who showed that viscosity can change the instabilities 
in the limit of very rapid conduction/viscosity, suppressing 
the growth rates of some of the modes (Fig. [TJ. 

In 2D simulations of the HBI, viscosity changes the non- 
linear saturation of the instability. Figure [3] shows that strong 
viscosity drives the magnetic field to bunch up in prominent 
vertical structures. However, in 3D the magnetic field lines 
slip past each other in interchange- like motions. These inter- 
change motions are not damped by parallel viscosity. As a 
result, in 3D the saturated state of the HBI (as measured 
by statistical quantities like (|6 Z |)) is not as strongly altered 
by the addition of viscosity, although it can take longer to 
reach the saturated state because viscosity slows the growth 
of the modes (Figure [¥J|. This conclusion is true even for very 
rapid conduction (and thus very rapid viscous damping), with 
Wcond(i) ~ 50ct>buoy (where the conduction frequency here is 
defined across the scale of the box L); this is a value charac- 
teristic of cluster cores (Fig. 0. 

In the cores of galaxy clusters the cooling times are short 
compared to the age of the universe. There is a long-standing 
problem of finding processes that can provide sufficient heat- 
ing to counteract this cooling. One possible solution is for 



conduction to bring in heat from the large thermal reservoir 
beyond the core. However, in the absence of other sources of 
turbulence, the HBI exacerbates the cooling flow problem by 
reorienting magnetic field lines to be more azimuthal, thus 
reducing the effective radial thermal conductivity, /s p . Our 
global simulations of realistic cool-core cluster models (Fig. 
[2} demonstrate that this reorientation of the magnetic geom- 
etry by the HBI is only minimally influenced by anisotropic 
viscosity (Fig. H & [TTJ . 

We also simulated galaxy cluster cores with additional 
turbulent forcing, intended to mimic the turbulence driven by 
galaxy wakes and/or AGN feedback. Turbulence is capable of 
suppressing the magnetic field reorientation of the HBI when 
the eddy turnover tim e is similar to the ins tability growth 
time (see Figure 11 of iMcCourt et al.ll2011al ). The interplay 
between turbulence, the HBI, and cooling in the cluster core, 
results in a very strong bimodality in the temperature pro- 
file and other cluste r properties such as the central entropy 
(|Parrish et al.l l201(t) . Modest levels of subsonic turbulence, 
less than 100 km s _1 , are capable of transforming a cool core 
cluster into a non-cool-core cluster. This turbulence is ener- 
getically small compared to the cooling luminosity in the clus- 
ter core; instead, the turbulence can be considered a catalyst 
that enables conduction to efficiently couple to the central 
regions. Larger turbulent velocities result in larger central 
temperatures and entropies, aided by both conduction and 
an inwardly-directed turbulent heat transport. With Bragin- 
skii viscosity, these conclusions about cluster bimodality are 
unchanged. In fact, we believe that viscosity aids the genera- 
tion of a bimodal cluster population: as the cluster core heats 
up, it becomes more viscous, which makes it somewhat more 
difficult for the HBI to re-orient the magnetic field, further 
promoting the transition to a non-cool-core cluster. 

Our numerical experiments are encouragingly consis- 
tent with the ob s erved bimodality in cl uster core properties 
(|Voit et al.ll200Sl ; ICavagnolo et al.ll2009r i. We find that cool- 
core clusters (with central entropies Ko < 30 keV cm 2 ) cannot 
be efficiently heated by conduction and must have an addi- 
tional source of heating. This source is plausibly AGN feed- 
back through bubbles, jets, or other coupling mechanisms. 
Indeed, cool-core clusters generically show evidence of ra- 
dio emission consistent with AGN activity. With turbulence 
above a critical threshold value, however, we find that con- 
duction is able to transform a cool-core cluster into a non- 
cool-core cluster with a much higher central entropy — this 
does not necessarily require AGN activity though it does re- 
quire roughly volume-filling turbulence. Observationally, this 
is consistent with the lack of significant AGN feedback indi- 
cators for high central entropy clusters. 

In local Cartesian simulations of the MTI we find that 
the evolution of kinetic energy is statistically unchanged by 
anisotropic viscosity; however, the amplification of the mag- 
netic field is diminished with the addition of viscosity (Fig. 
[6]). Magnetic fields are amplified by turbulence increasing the 
length of magnetic field lines; this motion is precisely the one 
damped by anisotropic viscosity. As a result, the magnetic 
field is amplified more slowly than the kinetic energy with 
anisotropic viscosity. For simulations with initially horizontal 
fields, increased viscosity leads to a small vertical bias in the 
magnetic field direction relative to isotropy. Simulations that 
start with initially vertically magnetic fields, a nonlinearly 
unstable configuration, evolve much more slowly at high vis- 
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cosity (bottom panel of Fig. [6]). This mechanism explains the 
small vertical bias seen for higher viscosity. 

In the outskirts of galaxy clusters, the temperature de- 
clines with increasing radius and the plasma is unstable to 
the MTI rather than the HBI. At these radii, the MTI can 
generate significant turbulent velocities up to an rms Mach 
number of (M) ~ 0.35 near r2oo- This turbulence represents 
a non-thermal source of pressure that can bias hydrostatic 
estimates of cluster masses. As emphasized in Parrish et al. 
(2011), this turbulence is in addition to that generated by 
structure formation. In the global cluster models of the MTI 
presented here (Fig. [HI), we find that viscosity only weakly 
affects this turbulence, as predicted analytically by Kll. The 
rms velocities with viscosity are only 5% smaller than the in- 
viscid case. Consistent with our local simulations, we find that 
the amplification of the magnetic field is modestly reduced for 
simulations with viscosity. The fact that anisotropic viscos- 
ity specifically damps all sources of turbulence that couple 
to magnetic field amplification makes it more difficult to in- 
voke turbulent mechanisms for amplifying the magnetic field 
in clusters from primordial values. 

Although our results with anisotropic viscosity do not 
significantly change any of the previous conclusions about 
the role of the HBI and MTI in the ICM, we believe that 
it is important to have established this conclusively. Most 
importantly, a fluid with a low collisionality like the ICM or 
the solar wind has a high viscosity, as the viscosity scales 
inversely with the collision frequency; the Reynolds number 
is correspondingly rather small (Fig. [2]). Thus the importance 
of viscosity is not a priori clear and needed to be explicitly 
studied. In future work it will be interesting to consider the 
role of anisotropic viscosity on other problems in the ICM, 
such as cold fronts or the large surface brightness fluctuations 
(ripples) observed in Perseus llFabian et al.| [2006). 



APPENDIX A: ALGORITHM AND TESTS OF 
ANISOTROPIC VISCOSITY 

Al Method 

We now describe the details of our algorithm for anisotropic 
m omentum transport. T his algorithm is similar to that used 
m iDong fc Stonel |2009l ). For simplicity consider the momen- 
tum and energy equations with only anisotropic viscosity 



d(pv) 
dt 



+ v-n = 0, 



0, 



(Al) 
(A2) 



where the viscous stress tensor is given by Equation [9] 
We proceed to discretize this in a very similar man- 
ner to how anisotropi c conduction was first discretized in 
iParrish fc Stonel |2005l ). For example, the first component in 
square brackets of the tensor has terms that can be written 
as 



bb : Vv 



f dv x £ dv x - dv x 
bx h by — h b z — — 

ox dy dz 



+ b y (...) + b z (...), 
(A3) 

where the ellipsis indicates similar derivatives of v y and v z , 
respectively. We will see momentarily that this stress is com- 
puted as an area-averaged, face-centered quantity which then 
requires velocity derivatives to be appropriately centered. 
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Figure Al. This figure illustrates the centering of quantities nec- 
essary for calculating the the shear on the x-face of cell (i,j) due to 
the transverse velocity gradient, dv x /dy. The velocity derivatives 
at cell corners (red arrows) are simple averages in Equation IA5I 
These derivatives are then averaged with a limiter (C) using Equa- 
tion lA4l to get the face-centered shear (thick black arrows) used for 
calculating the stress. 



The calculation of b x dv x /dx is straightforwardly dis- 
cretized from the cell-centered velocities. The calculation of 
the transverse velocity gradients, such as b x b y dvx/dy requires 
special care as monotonicity is required. This issue surfaces 
in anisotropic thermal conduction when it proves essential 
to limit the slope of the transverse temperature derivatives, 
e.g. the term b x b y dT/dy in the calculation of q x , to pre- 
vent heat from being con ducted from cold to hot regions 
(Sh arma fc Hammettl T2007'). This type of non- monotonicity 
can lead to negative internal energies with conduction. The 
same symptom can occur with viscosity, and it is necessary 
to ensure that momentum flows in the correct direction at all 
times. We therefore interpolate the transverse velocity gradi- 
ents to cell faces with a slope limiter: 



dv x 
dy 



C 



dv x 
dy 



dvx 
dy 



(A4) 



where i — i represents the left face of cell and C repre- 

sents a slope limiter. The velocity gradients at cell edges are 
simply arithmetic averages, e.g. 



dvx 
dy 



x,i — l,j-\-l 



Ay 



! !i -| 



Ay 



(A5) 

We choose the monotoni zed central dif ference (MC) limiter 
for our slope limiter (see Ivan Leerll 19791 . for details). 

Our explicit methods for anisotropic conduction and vis- 
cosity both have timesteps that are more restrictive than the 
MHD timestep as they are parabolic equations. The momen- 
tum diffusion has a Courant-limited timestep that is propor- 
tional to (Ax) 2 . Therefore, we choose to sub-cycle both dif- 
fusive operators with respect to the MHD timestep. For the 
anisotropic viscosity, we are able to simply call the standard 
MHD boundary conditions in between each sub-cycle. Due 
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to the smaller physical diffusivity, we are able to take much 
longer timesteps for viscosity than conduction. Sub-cycling 
the conduction and viscosity do not change the evolution of 
the HBI compared to reducing the MHD timestep to match 
the conduction timestep. 

After these preliminaries, we are now prepared to out- 
line the method for anisotropic momentum transport. Careful 
attention is paid to the centering of the terms in the differ- 
ence equations. Our method resembles the underlying Go- 
dunov scheme for Athena in that conservative cell-centered 
variables, e.g. momentum, are updated by differencing face- 
centered fluxes, e.g. the viscous stress — such a scheme is man- 
ifestly conservative. Our method for calculating one sub-cycle 
of a cell centered on is as follows: 

(i) Apply the MHD boundary conditions to synchronize all 
ghost zones across processors. 

(ii) Calculate face-centered velocity gradients, including 
limiters on transverse gradients, e.g. (dv x /dy) i _i ^ 

(iii) Calculate the face-centered viscous prefactor 

(Hi-ij = \ \[PV{T, P \_ X . + [pv(T,p)] i+lij } , (A6) 

where the diffusion coefficient (sometimes known as the kine- 
matic viscosity) depends on density and temperature as v oc 
y5/2^-i rp^g density dependence drops out of the stress ten- 
sor. The density dependence, however, is relevant for calcula- 
tion of the viscous timestep. 

(iv) Calculate the face-centered viscous stresses, e.g. 

n i on the left s-face. 

1 2 •■> 

(v) Difference the viscous stresses to get a cell-centered 
update to the energy and momentum equation, which is in 
the rc-direction: 

A(p»«) (ii = ^ (n i+ i/2,j - n 4 _ x/a>i ) , (A7) 

A(£) - = ^[( n -«W 2 „-( n '<-i/ 2 J- (AS) 

There are similar terms in the y- and z-directions. The At 
here is the sub-cycle timestep. 

Steps (i)-(v) are repeated for each sub-cycle. 
A2 Verification 

It is important to verify a new numerical method, such as 
anisotropic viscosity, with test problems with known solu- 
tions. To wit, we begin with two tests of the physics of linear 
MHD waves. We take advantage of the property that Bra- 
ginskii viscosity damps viscous motions only along field lines; 
however, our numerical implementation will have some spu- 
rious u±. 

We test our algorithm by measuring the viscous damp- 
ing rates of MHD waves. These damping rates are sensitive to 
both the magnitude and anisotropy of the viscosity and there- 
fore provide a stringent test of our code. Moreover, when the 
viscosity is small enough (quantified later), the damping rates 
can be calculated analytically, permitting a rigorous point of 
comparison for our simulations. In order to test the multi- 
dimensional nature of the code, we pick the wave vector k 
to be 45° from the grid axes and the magnetic field B to be 
56.3° from k. Thus, nothing is aligned with the grid axes and 
all terms in our algorithm are evaluated. 



Alfven waves directly probe the anisotropy of our algo- 
rithm. The fluid motions induced by these waves are purely 
transverse to the magnetic field and therefore are not damped 
by parallel viscosity. Any measured damping represents a spu- 
rious perpendicular (or isotropic) viscosity due to errors in the 
algorithm. We confirm that our code does not damp Alfven 
waves: even with an extremely large viscosity v — 1A 2 /T, 
where A is the wavelength and T is the wave period, we mea- 
sure a damping rate of 1.8 x 10 _4 /T. This damping is only 
~ 40% larger than that produced by the inviscid, default ver- 
sion of Athena with no explicit viscosity. For more reasonable 
viscosities, this damping is much reduced. Thus, the spurious, 
cross-field diffusion in our algorithm is typically much smaller 
than other sources of numerical dissipation in our simulations. 

We test the parallel component of the viscosity by mea- 
suring the damping of fast and slow magnetosonic waves. 
Both of these waves have a component of velocity parallel 
to the magnetic field and thus damp at a rate 

a = -^[(k-v)-3(b-k)(b-v)^ , (A9) 

where b, k, and v are unit vectors in the directions of the 
magnetic field, wave vector, and wave perturbation veloc- 
ity, respectively (|Braginskiilll965l ). Figure IA2I compares the 
damping rates measured in our simulations with this analytic 
expectation. In the left panel, we plot the ratio of the mea- 
sured to expected damping rates. The disagreement at low 
viscosities is most likely caused by the numerical, isotropic 
viscosity inherent in Athena's integrator. This viscosity is not 
included in the analytic expectation, so equation |X9] underes- 
timates the damping rate in this limit. With higher viscosity 
v > 0.1, the linear damping assumption in the derivation of 
equation I A9I is not valid. The algorithm performs well in be- 
tween these limits. 

In order to test the the numerical convergence of our algo- 
rithm, we chose a fiducial viscosity v = 0.01A 2 /T , such that 
the damping remains small but is not dominated by numeri- 
cal viscosity. We plot the fractional error in the damping rate 
Aa/a = (a mcas - <7 thcory ) / '(Theory in the right panel of Fig- 
ure lA2l This is not the usual L2 norm that is often plotted for 
linear wave convergence tests. Even at 16 zones/wavelength, 
the damping rate is correct to order unity. The quantity Aa/a 
converges at third order. This excellent convergence results 
from a convolution of the increasing accuracy of the linear 
waves themselves as well as the improved accuracy of the vis- 
cous transport, both of which are second-order algorithms. 

Perhaps the most demanding test of our combined al- 
gorithm for MHD, anisotropic conduction, and anisotropic 
viscosity is our measurement of the linear growth rates of the 
HBI and MTI in Figure [1] We are able to confirm the analytic 
growth rates to within ~ 1%, but unfortunately cannot reach 
the same accuracy as with the linear MHD waves. 

At least two effects limit our precision when measuring 
the growth rates of the HBI and MTI. First, the analytic 
growth rates to which we compare are only strictly known in 
the WKB limit kH 3> 1. Practical considerations limit our 
simulation setups to kH ~ 100, however. We expect correc- 
tions to the analytic growth rates of order l/(kH) ~ 1%, 
similar the typical discrepancy in our measurements. Addi- 
tionally, the structure of the eigenmodes of the HBI and MTI 
depend on the growth rate. Since we only know the growth 
rate to an accuracy of ~ 1%, we cannot initialize the simula- 
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Figure A2. Measurement of the damping rates a of the fast and slow magnetosonic modes. (Left:) The ratio of measured to predicted 
damping rates as a function of viscosity (in units of wavelength 2 / wave period). The agreement is very good over ~ 2 orders of magnitude 
in v. The disagreement at low u is likely caused by numerical viscosity in other parts of Athena's algorithm, while the disagreement at 
large v represents a breakdown in our linear damping calculation of a. (Right:) The fractional error in the analytic damping rates for our 
fiducial viscosity, uT/\ 2 = 0.01, as a function of resolution. The grey line represents third-order convergence. This excellent convergence 
result is the convolution of the second order accuracy of the linear waves themselves with the second order accuracy of our viscous transport 
algorithm. 



tion in an exact eigenstate. Thus, the perturbations we apply 
do not grow strictly exponentially in our simulations; there 
is an initial period of order £buoy during which the perturba- 
tions settle into their respective eigenstates. In practice, this 
limits the period of growth from ~ 2£buoy , when exponential 
growth begins, to ~ 4fb UO y, when the instability begins to 
become nonlinear. Thus we have only a limited time window 
over which to fit the data and cannot reach an arbitrary level 
of precision. It might seem that we could delay saturation, 
and therefore improve the accuracy of our growth rate mea- 
surement, by using a smaller initial perturbation. In practice, 
however, our boundary conditions do not hold hydrostatic 
equilibrium perfectly and eventually generate motions of or- 
der 5 x 10 _5 c s . Thus, when we use a weaker perturbation we 
lose it in the noise. 

Despite these limitations, we are able to confirm the 
growth rates with a comparable accuracy to the validity of 
the analytic theory. This, along with our measurements of 
the damping of linear MHD waves, represents a strong vali- 
dation of our code. 
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